#ifndef ORTHANTSAMPLER2_HPP
#define ORTHANTSAMPLER2_HPP
//#include "Data/Data.hpp"
//#include "Gibbs.hpp"
template<class Kernel, class Resample>
class Sampler2{
public:
	Sampler2(Kernel *K,Resample *R,Orthant2::Orthant *D,int M, double C){
		_p=(*D).Get_p();
		mat Xt(M,_p);//p le nombre de parametre a obtenir de Density
		W=new double[M];
		X=Xt;
		_K=K;
		_R=R;
		_Z=0;
		_C=C;
		_boo=0;
		_D=D;
		mat m(_p,1);
		m.fill(0);
		mat Sig(_p,_p);
		Sig.fill(1);
		Sig=diagmat(Sig);
		_G=new Distribution::Gaussian(_p,m,Sig);
		kk=-1;
		_thres=C*M;
		cout <<"thre " << _thres;
		_M=M;
		Normalize();
		y=X.t();
	//	cout << "///" << _n << "///";
	}
	~Sampler2(){
	//	delete[] W;
	}
	void Filter()
	{
		this->Init();
		int i=0;
		double p=0; 
		//for(int i=1;i<(_n);i++)
	//	i{
		y=X.t();
		_D->Set_i(i);
		_D->Weight(y,W);		
		Ess=Ess_W(W,0);
		cout << " ESS: " << Ess;
		X=y.t();
		for(int i=0;i<_p-1;i++){
	//		cout << i;
			//cout << " step: " << i << "\n";
			this->Step(i);
		}
		Ess=Ess_W(W,1);
	//	Correction2();
		
	}
	void Init()
	{
		//pior
		X=_G->scrambled(_M);
	}
	void Step(int i)
	{
		cout << "Ess:" << Ess << "\n";
				//cout << "\\";
		if(Ess< _thres)
		{
			double ess=Ess_W(W,1);
			mat v=growingvect(_M);
			double *w=new double[_M];
			for(int j=0;j<_M;j++)
			{
				w[j]=exp(W[j]);
		//		cout << w[j] << " ";
			}
		//	cout << "\n";
			
			(*_R)(&_M,w,v);
			delete[] w;
			//cout << v << "\n";
			Arangemat(X,v);
			Normalize();//poids fixé a 1
			//mat C=cov(X);
		//	Normalize2();//poids fixé a 1
		//	for(int j=0;j<_M;j++)
		//	{
		//		cout <<exp(W[j]);
		//	}
			
			//cout << C;
			//cout << "sum: "  << sum(X,0) << "\n";
			//cout << "before:" << X.row(1);
	//		for(int k=0;k<3;k++)
	//		{	
			
			//cout << "\nmean " << mean(X,0) << "\n";
		//	cout << "\nmean " << mean(X,0) << "\n";
/////////////////////////////////////////////////////////////////////////////////////

			for(int j=0;j<1;j++)
			{
				_K->Move2(&(X),i);
			}
		}
		//	cout << "\nmean " << mean(X,0) << "\n";
		
			y=X.t();
			//cout << X << "////";	
	//////////////////////////////////////////////////////////////////////////////
		//	b+=(double)(1/(double)(_P-1));
			_D->Set_i(i+1);

			//cout<< "before:" << W[0] ;
			//cout<< "after:" << W[0] ;
			//dans le cas de tempering peut etre determiner
			_D->Weight(y,W);//le p designera comment on avance 
		
			
		/*	for(int i=0;i<_M;i++)
			{
				cout << W[i];
			}
			*/
			//	cout << "after:" << y.col(1);
			X=y.t();
			Ess=Ess_W(W,0);
			//cout << " ESS: " << Ess;
							
			//cout << " ESS: " << Ess_W(W,1);
	}
	double Ess_W(double *w, int tf)//renormalisation???
	{
		double *foo2= new double[_M];
		double sum2=exp(w[0]);
		for(int i=1;i<_M;i++)
		{
			sum2+=exp(w[i]);
			//if(tf)
			//{
				cout << exp(w[i]) << " ";
			//}
		}
		cout << "//////////////////////////////\n";

		if(tf){
			_Z+=log(sum2)-log(_M);
			cout << "Z "<< _Z << "\n";
		}
		for(int j=1;j<_M;j++)
		{
			foo2[j]=exp(w[j])/sum2;
		}
		sum2=foo2[0]*foo2[0];
		for(int i=1;i<_M;i++)
		{
			sum2+=foo2[i]*foo2[i];
		}
		delete[] foo2;
		
		
		return (double)1/sum2;

	}
	
	mat Get_theta(void){return X;}
	void Normalize(void){
		for(int i=0;i<_M;i++){
			W[i]=0;
		}
	}
	void Normalize2(void){
		double sum=0;
		for(int i=0;i<_M;i++){
			sum+=exp(W[i]);	
		}
		for(int i=0;i<_M;i++){
			W[i]=W[i]-log(sum);
		
		}
	}
	mat   Get_W(void){
		double sum=W[0];
		mat C(_M,1);
		for(int i=1;i<_M;i++)
		{
			double t=log_add(sum,W[i]);
			sum=t;
		}	
		for(int i=0;i<_M;i++)
		{
			C(i,0)=exp((W[i]-sum)+log(_M));
		}
		return C;
	}
	double Get_Z(void){
		return _Z;}
	mat Get_Phiv(void)
	{
		int n=Phiv.size();
		mat res(n,1);
		for(int i=0;i<n;i++)
		{
			res(i,0)=Phiv[i];
		}
		return res; 
	}
	mat Get_Acc(void)
	{
		int n=Acc.size();
		mat res(n,1);
		for(int i=0;i<n;i++)
		{
			res(i,0)=Acc[i];
		}
		return res; 
	}

private:
	mat X;//value particle
	double *W;//weigths
	double Ess;
	int _M;
	double _thres;
	double _Z;
	Kernel *_K;
	Resample *_R;
	Orthant2::Orthant *_D;	
	mat y;
	double _temp;
	int _p;
	int kk;
	vector<double> Phiv;
	vector<double> Acc;
	Distribution::Gaussian *_G;
	int _boo;
	double _C;
};



#endif
